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Abstract. We present a theory of frequency-dependent counting statistics of 
electron transport through nanostructures within the framework of Markovian 
quantum master equations. Our method allows the calculation of finite-frequency 
current cumulants of arbitrary order, as we explicitly show for the second- 
and third-order cumulants. Our formulae generalize previous zero-frequency 
expressions in the literature and can be viewed as an extension of MacDonald's 
formula beyond shot noise. When combined with an appropriate treatment of 
tunneling, using, e. g. Liouvillian perturbation theory in Laplace space, our 
method can deal with arbitrary bias voltages and frequencies, as we illustrate 
with the paradigmatic example of transport through a single resonant level model. 
We discuss various interesting limits, including the recovery of the fluctuation- 
dissipation theorem near linear response, as well as some drawbacks inherent of 
the Markovian description arising from the neglect of quantum fluctuations. 
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1. Introduction. 

Transport of electrons through nanoscopic conductors is a very powerful tool to learn 
about interactions and to characterize quantum systems [1]. Examples include the 
quantum Hall effect [5], weak localization [3], or universal conductance fluctuations 
[4]. Transport processes are governed by tunneling events, which are stochastic in 
nature. It is therefore natural to expect that the statistics of these tunneling events 
will be strongly influenced by interactions and quantum effects. Interestingly, these 
statistics, which can be analyzed by studying current fluctuations, contain a great 
deal of new information beyond that provided by dc transport [SJ [S]- In particular, 
the second-order current correlation function (noise), can be used to determine the 
effective charge [7] and the statistics of the quasiparticles [8] , and to reveal information 
on the transmission properties of the conductor [5] . Moreover, current correlations can 
be used to learn about entanglement [9l [10] , quantum coherence [11] , and the deep 
connection that exists with fluctuation theorems [121 [El [H]- Further information 
can be gained from noise at high frequencies which is valuable for extracting internal 
energy scales in systems such as quantum dots (TS] , spin- valves [16] , Cooper Pair Boxes 
[17], diffusive wires [18] or chaotic cavities [T9] . 

While most of the work on noise is theoretical (in particular at high frequencies) , 
the field of noise and counting statistics is producing a great deal of experimental 
breakthroughs, including measurements of high order cumulants I^Tl \n[ 
[23 [ini [371 [2H1 [IS]- Owing to this experimental progress, noise measurements at high 
frequencies, which until recently were scarce, are now possible [SO] [311 [32l [33l [34l l35l 
[361 [371 [38]. 

A proper treatment of fluctuations in non-equilibrium transport is needed to 
address the problems listed above. While the list of theoretical available is too large 
to be given here, it is safe to say that they can be divided roughly into three families: 
the scattering approach [39] [5] , the Keldysh Green's functions method [40] , and the 
various quantum master equation (QME) treatments |41| (for a recent overview of 
transport in this last context see e.g. Ref [42]). A theory of counting statistics of 
electron transport was first formulated by Levitov and Lesovik for noninteracting 
electrons using the scattering formalism |43j . and later works enabled the treatment 
interacting problems [44] [45] . QME approaches followed [46] , and proovcd particularly 
useful for studying systems in the Coulomb Blockade regime. Recent avances within 
this last scheme also involve studies of the counting statistics including non-Markovian 
dynamics [47l [48l [49l [50] . 

In Ref. [HI] we developed a method for calculating high-order current correlations 
at finite frequencies in the context of Markovian QMEs, and the aim of this paper 
is to significantly extend that work. We provide a detailed derivation of our multi- 
time generating function, Eq. ([27| . We present a new approach to derive finite- 
frequency cumulants from this expression, noticing that only the stationary solution 
of the problem is required. We give analytical expressions for the second and third- 
order current cumulants (noise and skewness respectively). These results generalize 
previous zero-frequency expressions in the literature and recover the finite frequency 
shot noise expressions [52l [53] obtained using the MacDonald's formula [54]. Our 
method can thus be viewed as a generalization of this formula, as it allows us to obtain 
high-order spectra such as the frequency-dependent skewness (Eq. ([35])). To illustrate 
the formalism we study the case of a single resonant level (SRL) model, and compare 
it with the exact solution and with the non-equilibrium version of the fiuctuation- 
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dissipation theorem, derived in various approaches, such as for tunnel junctions, or 
for the weak cotunnehng regime in quantum dots [551 IMl EZl \SS\ ■ 

The paper is organized as fohows. In section [5J we present our formahsm of 
finite-frequency cumulants in the context of Markovian QMEs. Subsections 12.11 and 
12.21 are devoted to establish the general framework of full counting statistics. In 
subsection 12.31 we derive a multitime cumulant generating function. Subsection 12.41 
shows how to obtain finite-frequency cumulants of the current distribution. Here, 
exact equations for the frequency-dependent noise and skewness are given. We end 
this part with special emphasis on how to calculate the counting statistics of the 
"total" and "accumulated" currents (subsection 12. 5p . We explicitly show how both 
current correlations and charge correlations can be calculated. In section [3l we study 
the example of a SRL model, providing spectra for the frequency-dependent noise 
and skewness and a detailed comparison with the fluctuation-dissipation theorem, 
the finite-frequency version of the non-equilibrium fluctuation-dissipation theorem, 
and the exact solution of the SRL model. First we focus on the zero-frequency case 
(section |3.1|) . where the general behaviour of noise and skewness is presented as a 
function of different system parameters. In section [3?2l we extend this study to the 
finite-frequency case. Interestingly, we show that, even though the theory does not 
contain quantum fluctuations, the Markovian limit is basically exact in transport 
configurations, level within the bias voltage window, as long as ftw ^ eV or huj ^ eV. 
In intermediate situations, where Huj ~ eV, or with the level outside the bias windows, 
the Markovian limit fails at finite frequencies due to its lack of quantum fluctuations. 
We also demonstrate that the noise spectra for particle currents and the ones for total 
currents significantly deviate from each other, even for large asymmetric coupling to 
the leads, namely Th/Tl ^ 1- In section U] we summarize our results. Most of the 
technical details and intermediate steps of the derivations in section [2] are discussed 
in detail in [Appendix A[ where we also present a diagrammatic technique to arrive 
to the expressions for the cumulants shown in sec. 12.41 In section [Appendix B| we 
describe how to calculate the kernel of the QME to lowest (sequential) order using 
perturbation theory in the Liouville space. 



2. Theory 

2.1. Quantum master equation 

We are interested in phenomena which can be described by an equation of the type 

m = ^p{t), (1) 

where C is the so called Liouvillian, governing the evolution, and p is the density 
operator of the total system. Specifically, our theory will be useful to processes 
amenable to the counting of a classical stochastic variable n, which can be, for example, 
the number of particles that have undergone a particular process in the system. We will 
focus in particular in transport systems, consisting on a central region, with a known 
set of many-body eigenstates {|a)}, and respective eigenenergies {Ea}, attached to 
non-interacting electronic leads at different chemical potentials. This set-up can be 
described by a Hamiltonian that takes the form H = Hs + 'Hr + Ht, where Hs and 
Hr refer to system and leads respectively, and Ht is the coupling between them. The 
different terms can be written as 

^5 = 5]i?a|a>(a| (2) 



Finite-frequency counting statistics of electron transport: Markovian Theory 



4 



■Hii = ^ (e,,Q + /ia) cl^Crja (3) 
Ht = ^ Vrianicljadm + H.C. (4) 

Here c]^^ creates an electron with quantum numbers ry in lead a, and dm annihilates 
an electron from site m in the central region, e^a are the eigenenergies of the electrons 
in the lead a, Vrjam is the coupling energy between a state in contact a and the level 
m in the system. the chemical potential of lead a, and that allows the system to 
be driven out of equilibrium. 

Given the Hamiltonian, the full evolution can be written in terms of the Liouvillian 
by using the von Neumann equation 

m = -^[n,p{t)]^cp{t). (5) 

We are actually interested in the dynamics concerning the central system. We therefore 
trace out the reservoir degrees of freedom, arriving at an equation of motion for the 
reduced density operator ps{t) = Trp{{p(i)}, that, if we assume the dynamics of the 
reservoirs to be much faster than that of the central system, can be approximated as 
a Markovian master equation 

Ps{t)^Wps{t), (6) 

where W is a kernel driving the dynamics of the system. The charge flow through the 
conductor is governed by the stochastic hopping of electrons in and out of the central 
region. This processes are susceptible to classical counting and thus the reduced 
density operator can be unravelled into components psin^, t), corresponding to having 
Ha = ...,—1,0,1,... extra electrons in lead a [SHI ED] ■ The kernel can also be split as 
W = Wo + S± where refers to the physical process in which one electron in 
created (+) or annihilated (-) at lead a, and Wo corresponds to the part in which no 
tunneling processes take place in that lead. It can be shown that ps [ua , t) fulfills the 
equation (see e.g. [6T]): 

PsK,t) = Wops(n„,0 + W+PsK - l,t) + W^psK + (7) 

valid provided that only single particle tunneling processes occur. Although Eq. ([7]) 
focuses on the single counting at a particular lead a, it can be generalized to account 
for tunneling processes of k particles at the different system-reservoir junctions: 

ps{ni, . . .,nM,t) = WoPsini, . . .,nM,t) 

+ ^ W^'''ps{ni,...,naTk,...,nM,t) (8) 

where Wq is the part in which the number of particles is not changed in the central 
region, M the number of leads, and k labels the process in which k particles "jump" 
at a time. 

Unfortunately, solving equation (|5]) in the n-space requires truncation to a certain 
n and diagonalization of a tridiagonal matrix. It is therefore more convenient to solve 
it taking its Fourier transform. Multiplying (|8]) by e™^-^^ . . . e^"'^'^'^' and summing 
over ni , . . . , we obtain 

Psix,t)=yV{x)psix,t), (9) 
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Here the counting field x refers implicitly to all counting fields, and we have 
defined ps(x,i) = E„i,...,„„ e^^x^ . . . e™^^»VsK, ■ • ■ , ^m, i) and W(x) ^ Wo + 
E^ _i_ yV^e^'-^°. For a time-independent kernel, the solution to Eq. (|9]) is 

ps{x,t) = n{x,t-tn)ps{x,to), (10) 
with the time evolution operator ri(x,^ — ^o) '■— e^'^-^'^*"*"-'. 

2.2. Full Counting Statistics 

Importantly, the knowledge of the system's density operator resolved in n allows us to 
obtain the the full counting statistics (FCS) of the system, that is the probability 
distribution P{ni, . . . ,nM,t) of the number of electrons transmitted through the 
system-lead junctions. This is accomplished by noting that 

P(7ii, . . . ,nA/,<) = TT{ps{ni,...,nM,t)}. (11) 

Transforming the probability distribution to the X"Space, we have the moment 
generating function (MGF) 

g{x,t)=Y,P{n,t)e"^^, (12) 

n 

where now n refers implicitly to ni, . . . , um- Using Eq. (fTQ|) one has 

gix,t) = TT{nix,t-to)ps{to)}, (13) 

equation that was already derived by Bagrets and Nazarov [35] . The iV-th. derivative 
of the MGF with respect to x gives the A^-th moment of the distribution of the number 
of particles that have gone in or out a particular lead a: 

.Nu.^^ 9^Six,t ), 

d{iXaY 

When equation ([T51) is used, averages with respect to the stationary state are 
established by taking ps{t{)) = pg*°* (defined by Wpg*"* = 0). This means that 
counting will start at a time in which the system has reached its steady state, 
and therefore the fluctuations we study are around this state Q. 

The moments of the current distribution can be calculated as [J 

{I^{t))^j^{<{t))- (15) 

This relation is important as it relates the stochastic variable n with the current 
of particles flowing through the system. Even though the current studied here is a 
classical variable, it contains quantum effects present in the system. In the formalism, 
these are inherited from the Liouvillian operator in equation ([l}. Generally, we are 
interested in the cumulants, rather than the moments, of the current distribution. 
These can be obtained from the derivatives of the cumulant generating function 
(CGF), defined by lnQ{x,t). Therefore we have 



{<it))^ J lx^o. (14) 



Ia)c — 37 a/- ^Tat' lx->0,t^oo, (16) 



d d^F{x,t 



f In the following, all the averages will be taken with respect to this steady state. An average in the 
Liouville space will be therefore written as (A) = ■ A • pg'"' = ((0|A|0)), where |0)) = p|'"' is the 
normalized stationary system density matrix (written as a vector), and tx = ((0| is the transposed 
trace vector that sums over the population degrees of freedom. 

I Throughout the paper we will use e (electron charge) = k (Boltzman's constant) = h (Planck's 
constant /27r) = 1. 



Finite-frequency counting statistics of electron transport: Markovian Theory 



6 



where {■■■)c denotes cumulant average |62] and the limit i — oo ensures that the 
average is performed in the stationary state. Also, notice that the probability 
distribution itself can be obtained by inverse Fourier transform of the MGF. 

From the x-indepcndent kernel of the reduced QME (jB]) , the x-dependence leading 
to equation (jO)) can be actually introduced in a simpler way than resolving the density 
operator in n and taking the Fourier transform. As we describe in [Appendix B[ it is 
enough to include counting fields in the appropriate tunneling terms of the Kernel |63j , 
and this procedure is fully equivalent to solving a generalized Von Neumann equation: 

p{x,t) - ~^{n+{t)p{x,t) ~ p{x,m-it)), (17) 

in which the time evolution in the forward (+) and backward (-) Keldysh part of 
the real time axis is governed by different Hamiltonians |55| . specifically "Hy = 
12ri a m^vame'^^^^^clj^dm + H.c. Tracing out the reservoir degrees of freedom in 
equation ([T7|) . one can get equation ([9|) and proceed to obtain the FCS of the system. 

2. 3. Finite-frequency full counting statistics 

In this paper we want to study correlations at finite frequencies, for which the 
scheme presented above has to be generalized. To this end, one needs to consider 
a joint probability distribution, P(ni, ti; . . . ; tiat, ^at) defined as the probability that 
ni electrons have undergone a particular process after a time ti, n2 electrons after a 
time t2, etc. Here we focus for simplicity on a particular lead, and denote n as the 
number of particles transferred to (from) it, with associated counting field x (~x)- 
Being straightforward to include processes at different leads. 

The connection between this joint probability and the density operator (analog to 
Eq. (Ilip ) is not straightforward. To connect them we first need to specify a prescription 
for the symmetrization of the cumulants and the probability distribution. This 
prescription actually depends on the detection scheme [Ml ESI ES] ■ Here we assume 
that "classical" detection is carried out, so the detector is incapable of distinguishing 
emission from absorption. This means that the results we will present correspond to 
the fully-symmetrized version of the power spectrum 

/oo 
dti... dtNe-"^''' . . . e-'""*" 
-OO 

xTs{I{tl)...I{tN))c, (18) 

where 7s is the symmetrization operator, that sums over all possible time (or 
frequency) switchings, that is, we have for example Ts{I{ti)I{t2)) = {d{t\)I{t-2)) + 

U(i2)/(il)). 

The spectrum psp can be derived from a iV-time (symmetrized) CGF 
defined by 

n\ ,njv 

where % := (xi, . . . , XAr)T, t [ti, . . . ,tN)T, n := (rii, . . . , n7v)T, P''^^ refers to the 
symmetrized joint probability, tj'^^ to the multitime MGF, and the subscript T to 
the transpose of a column vector. That is, we have 

/CJO 
dti... dtNe-'"'*' . . . e"''""*" 
-OO 

xai,...9t„a,^,...a,^^J-W[X,t] I^^Q. (19) 
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/>oo 

S^'^'>{uji,...,ujn) {iuji)...{iujn) / dti . . . dtAfe-'"^*^ . . . e"*'^"*" 



xa,^,...%„^W[x,<] l^^o- (20) 

Both the probabihty P(^) [n, t] and the CGF T^^'> [x, t] can be calculated from the 
density operator and the Kernel W if we use the Markovian approximation in the 
coupling with central system-reservoir coupling. Within that limit we have the 
evolution local in time given by (|10p and also the factorization property 

P>{ni,ti; . . .;nN,tN) = P{ni,ti)P{n2,t2\ni,ti) 

X . . . P{nN,tN\nN^i,tN-i), (21) 

where the symbol > constraints the times to tk > tk-i- Notice that as we 
are considering the totally symmetric correlation function, we need to consider 
P(^)[n,t] = TP^ {ni,ti; . . . ]nN,tM), where T is the time-ordering operator [5T] . 
P(n, t\n' , t') is the conditional probability of counting n electrons after time t provided 
that we counted n' electrons after time t' , and can be computed as 

?•'-';;;?'"'''» (22) 

Tripsin', t')\ 

where the normalization in the denominator accounts for the collapse of the state due 
to the measurement, as given by the von Neumann's projection postulate [67]. r2(n, t) 
is the propagator in the n-space, that is. 



PS 



{n, t) = Y, - t - t')Ps{n', t'), (23) 



and can be extracted from equation ([5]) or by inverse Fourier transform of the 
propagator in the x-space: 

n{n,t)^ J ^e-^-m{x,t) (24) 

An expression for the joint probability distribution in terms of propagators can be 
then derived using (PT|) together with and . Alternatively, it can be obtained 
using the Chapman-Kolmogorov property for Markovian evolutions, from which we 
have 

P{nN;tjy) =Tt ^ n{nN - nN-i;tN ~ tjy^i) 

Til ,...,TIJV-1 

X rt{nN-i ~ nN-2] tN-i - tN-2) ■ ■ • ^l{ni;ti)ps{to) (25) 
As we also have P(njv;tAr) = J2n n _ P"^ ("-i, ^i; ■ • ■ ; "-Wi ^Af); reminding that 
Ps(io) =pr*' wcfind 

JV 

P^^^[n,t]=r(l[n{,.N-k,rN-k)), (26) 
fc=i 

where Uk := rik+i — nk, Tk := tk+i~tk and (•) := Trj^pg*"*}. Transforming expression 
([26]) to the x-space, we find the CGF 

N 

J-W[x,<]=lnr(nf^fe,T^-fc)), (27) 
fe=i 
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being Xk X)i=Ar+i-fc Xi- The structure in Eq. (|27p is encountered in many branches 
of physics such as statistical physics and field theory, where connected correlation 
functions arc obtained by taking derivatives of the logarithm of the corresponding 
generating functional (the partition function, the S'-matrix, etc). Note in particular 
the analogy with the partition function presented in Ref. |62j . 



2.4- Finite-frequency cumulants 

Equation (|27p permits us to obtain frequency-dependent current cumulants to 
arbitrary order. This was used in Ref. [51] to study the second and third cumulant 
in single and double quantum dot models. Explicit derivatives of ()27|) and the 
eigen-decomposition of the Kernel was used then to that end. In this subsection 
we show that only the stationary solution of the problem (solution to an algebraic 
equation) is needed to compute the finite-frequency current cumulants. We give 
analytical expressions (valid within the Markovian approximation) for the noise 
(second cumulant) and skewness (third cumulant) of the distribution of charge fiowing 
through a conductor. 

Let us decompose the Fourier transform in equation (jlSp into a set of Laplace 
transforms (defined as f{z) := dte~^* f{t)), and the cumulant averages in terms of 
moments (c.f. for example Eq. (2.8) in Ref [HS])- Doing this we find[^ 

5W>(zi) = 5W>(zi), (28) 
S('>iz,,z,)^Sl?>iz,,z,)~^:^^ (-1) (/)^ (29) 

5(3)>(^l,Z2,Z3)=5(?)>(zi,Z2,Z3) 

'-1 

21 



/)5(2)>(Z2,Z3) 
— ) (/)5(?)>(Z1,Z3) 
— ) (/)5(?)>(Z1,Z2) 



2(-V-V-)W'- (30) 



Z\ ) \Z2 J V 23 

The notation ">" denotes the unsymmetrized correlation function corresponding to 
the time ordering t^ > ... > t2 > ti. Symmetrization in the frequency space 
implies adding the part corresponding to negative z and summing over all the possible 
switchings of frequencies. The subscript m means moment. These can be obtained as 

>(zi, . . . , zat) = zi . . . ZNd,^, . . . d,^,g^^^> [x, z] (31) 
with z = (zi, . . . , zn)t and[jil 

N ^ 

g^''^>[x,z] = (l[n{xk,Sk)) , (32) 

fe=l 

§ Notice that eqs. 1281 1- 11301 1 can also be derived if the derivatives of the CGF are decomposed in terms 
of derivatives of MGFs. For example, for AT = 3 we have J^jjs = ^123 "^f '^M -S2^'Sl3' 5i2' + 
2ef> g'^'g*'^^ with /, := axi/lxi=o, fij ■■= 9xi9xj/lxi,Xj=o. fijk ■■= 9x,9xj9xfc/lx,,Xj,Xfc=o- 
II In the frequency domain, the prescription > can be taken similarly, that is zjv > . . . > 22 > ^ii 
and finally symmetrize the result. 
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The advantage of having moment averages is that we can use a diagrammatic 
technique (see [Appendix A] ) to easily obtain the correlation functions. Symmetrizing 
5'(JV)>j2;] and evaluating it at 2 = iu} (being u) = (wi, . . . ,lon)t) we find that S'(^^[u;] 
is proportional to 5{oji + ... + lun), as required by time-translational invariance. 
Defining the jump supcropcrators J^^ := [W(x) — W(x = 0)] and their derivatives 
Jq^^^ := d^J^l^^Oi we arrive to the following expressions for the current, noise and 
skewness of the current distribution (c.f. [Appendix A| for details): 

(33) 



ilstat — (i7o ') 



{Jo 
{Jo 



(1) 



''r!o(-zc.')Jo''^) + (Jo^''f^o(-*^)jJ'^) 

''r!o(-zt^)Jo^''r!o(-*'^')^J''> 

^'^17o(iw')J'o'^f^o(«c^)J"J'^) 
^'f7o(-za;)J(J^'no(za;' - iu;)Jo^^') 
^^17o(ic^')J'o^^^o(iw' - it^)Jo^^^) 

'^QoiiuJ - iu;')J^''^noiiuj)J^''^), (35) 

being Qo{z) := [z — W{x = 0)]^^. These equations generalize the zero-frequency 
results found in Ref. [5T] (c.f. their eqs Eqs. (7) and (8)) to finite frequencies 
(see [Appendix A| for the zero- frequency limit of ((33l) - ((35l) ). Results for higher-order 
cumulants can be similarly obtained. 

The relation between cumulants and moments can be formally expressed 
more generally at the level of the generating function. To do this one should 
follow the derivation by Kubo [S^, making use of the property {expi^^niXi)) = 
exp{{exp{^^niXi) — l)c} in our context, arriving to a similar result to (7.25) in 
Ref. [62]. This allows for the calculation of frequency-dependent cumulants of the 
current distribution up to any order, reproducing in particular the results presented 
above. If a diagrammatic expansion in the Liouvillian space is used [681 1691 I70| . 
cumulant averages become particularly useful since one can then keep only connected 
diagrams as those contributing to the average, in a similar way that this is done in 
quantum field theory. 



2.5. FCS of total and accumulated currents 

At finite frequencies, to have a theory consistent with current conservation it 
is essential to include the so called displacement currents (S]. This point is of 
vital importance to reproduce correctly the noise spectra measured experimentally. 
Although our discussion has focused, by construction, on particle currents so far, we 
show here how to include the effect of displacement currents in our formalism. 

Let us illustrate this point by considering a quantum well with two terminals (L 
and R) in contact with Fermi leads at different chemical potentials. There will be 
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then a net current flowing through both terminals, but also, charge can "accumulate" 
in the well for some time. Therefore charge conservation can be expressed as 

Q{t) ^ II — Ir = laccum, (36) 

with Q the charge in the well and II, Ir referring to the currents through the left 
and right contacts respectively. Q represents the displacement current, Idis^ which 
can be partitioned as Idis = {a + [i)Idis = I^is + ^dis^ where a and f3 describe how 
the displacement current is divided between the right and left reservoirs (obviously 
a + p = 1). This partitioning allows us to write the current conservation as 
II - li, - {Ir + I^i^ = 0. Equivalently, hot^ h - iL = Ir + = ^II + I3Ir, 
which is the so-called Ramo-Shockley theorem (jj. In the simplest wide-band limit, the 
partition coefficients can be written in terms of tunnel rates only as a := ^r/{^l+^r) 
and (3 := ^l/(^l + ^r) I7I]0, and this will be the partitioning we will use throughout 
the paper. 

Experimentally, one can measure correlations of the current through the device by 
transport measurements [20l[21], or indirectly by studying the current through a charge 
sensor, such as a quantum point contact [22l [23l [24], that reveals whether the well is 
"charged" or "uncharged" . The second method gives the statistics of the transport 
current only for very large bias voltages (unidirectional counting) but, in general the 
time-dependent transport current and the charge statistics are different. Morevover, 
when the device itself is used as a detector, the difference between transport 
fluctuations and charge fluctuations leads to profound physical consequences. Unlike 
the inelastic backaction induced by current fluctuations of the detector |65j , the one 
induced by charge fluctuations is the fundamental Heisenberg backaction associated 
with the measurement |73j . Both transport and charge fluctuations can be accounted 
for in our formalism by considering counting fields 

Xtot ■■= XL + XR, (37) 

Xaccum ■■= PXL " OiXR, (38) 

which lead to respective jump operators 

j(:i^a"^S+/^"^S' (39) 

jin) ^ 7'(") I 7'(") (AO) 

where J7J^,l and J^x,R refer to the two independent tunneling processes occurring at 
each barrier. 

The "total" cumulant through a two terminal device can be then calculated 
performing derivatives of the CGF with respect to Xtot defined in ([57)) . This leads to 
expressions ([55|) - ([55|) with J'q"^ substituted by JqII^. everywhere. Also, the spectrum 
of charge fluctuations 

/•oo 

slf\uJi,...,ujN) := / dti...dtNe-"^''' ...e-'^"*" 



follows from 



■^Ts{Qiti)...Q{tN))c, (41) 

(lUJi)...[lUJN) 



^ In this paper we mean "total" cumulant when a subscript is ommited. 

+ In a Coulomb Blockade model, the partition cocfBcients read a = Cr/{Cl + Cn) and /3 = 
Ci^KPij -l-Cfl), where C/,, Cjj are the capacitances of each barrier and we have neglected capacitive 
effects from the gate. See for instance | 72| . 
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with Saccum[<-^] obtained from cqs. upon the change Jq — !• J^o,acumm- For 

example, Sq\ll>) — (l/a;^)5acctim('^)- Notice that for a capacitive conductor, due to 
the relation between charge and vohage, this charge noise is proportional to the voltage 
noise. Finally the "left" and "right" cumulants can be computed with Jo — >■ J7o,l and 
Jo Jo,B. respectively in (|33 |) - (f35|) . 



3. Results 



To illustrate our method, we analyze the transport statistics of the prototypical 
example of spinless electrons passing through a SRL model. The system consists 
on a two-terminal conductor with a discrete energy level in the central region, and is 
described by the Hamiltonian 

n^e\l){l\+ (efea+Ma)cLcfea+ E V,„cL|0)(l| +H.C., 

k,a£L,R k,aGL,R 

(43) 

where k is the momentum and |0) and |1) are the only two possible states (referring to 
empty and occupied level) due to Coulomb blockade, with respective energies and 
e. 

In the infinite bias limit (voltage V much larger than the other energy scales, 
excepting the bandwidth of the Fermi leads) the Hamiltonian (|43p leads to the kernel 

expressed in the basis {|0), |1)}, and where Fq ~ Ta{E) := ^ J2k \^ka\'^S{E — eka) are 
the rates accounting for the system- reservoir coupling. Using and the 

simplicity of the model allows us to derive analytic results in this limit; for example 
the current gives Istat = ^l^r/^, where F := Tl+Th, and the "total" noise expressed 
in terms of the Fano factor (i^^^' := S'-^yigtat) reads 

ri+r| + (l-2a^)u;2 

At finite bias voltages, the kernel in Eq. is no longer valid. Among the 

various choices to calculate W(x) in this case, we use Schoeller's approach [68l [691170] 
(c.f. [Appendix B[ ) which allows us to calculate the Markovian kernel to the desired 
order (sequential tunneling in our case) without further uncontrolled approximations 
(such as the secular approximation). It is important to mention that the frequency- 
dependent shot noise of the SRL model can be solved exactly [71], and therefore one 
does not need to use the above approximations. However, to the best of our knowledge, 
a finite-frequency study for this model beyond shot noise is yet lacking. Here we use 
the exact solution as a benchmark of the Markovian approximation in order to identify 
the regions of validity of our theory. This benchmark is important because most of 
the papers in the literature discussing shot noise in the context of QMEs make use of 
the Markovian approximation. 

Another important check for the theory is to reproduce the fluctuation-dissipation 
theorem (FDT) in the appropriate regimes. Near linear response, that is, for applied 
voltages V much smaller than the temperature T, the low-frequency noise spectrum 
should follow the Jonhson-Nyquist relation S^^^ = 2kTG [75l [Tg, where G is the 
dc conductance. This equilibrium FDT was later extended by Callen and Welton to 



f(2)(,)^il±iiLtii_i^W^. (45) 
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include quantum fluctuations f77| . relevant when the measured frequencies are larger 
than the temperature. The FDT takes then the form S'^^^oj) = fujcoth{^)G{uj), 
where G(w) is the ac conductance. This expression can be equivalently written in 
terms of the Bose-Einstein distribution Nbe{^) = l/[e^ — 1], since coth(^y) = 
2Nqe{^^) + 1 = NbeI^) — NbeI^^^), and it becomes clear that the symmetrized 
noise, which we are considering here, contains both absorption and emission. Out 
of equilibrium, a fluctuation-dissipation relation can be also found for some particular 
cases, such as tunnel junctions |561 157| or for quantum dots in the weak cotunneling 
regime |58j . For a two-terminal conductor driven out of equilibrium, the symmetrized 
noise spectrum takes the general form [751 IZHl 1311 IS 

^ " 2 \ 2kT J 

^A.(l-i?„), (46) 

n 

where D„ is the transmission coefficient of the conduction channel n. This expression 
has various interesting limits. First, in the tunneling regime (D„ <C 1), it gives the 
non-equilibrium fluctuation-dissipation theorem (NEFDT) as reported in [SSI HZ] for 
tunnel junctions: 




expression that is also exact for quantum dots in the weak cotunneling regime |58j , and 
whose zero-frequency limit 5*'^' = IstatCoth{^^) has been derived in the context of 
counting statistics [55]. For low voltages, eV ^ kT, equation (|46)) recovers the Callen 
and Welton equilibrium relation, and if also huj ^ kT, it gives the Johnson-Nyquist 
FDT (thermal noise regime). Finally, if eV ^ kT, huj (shot noise regime), we find 
5*^^) = (■/, where the coefflcient ( = '^^^ Dn{l ~ Dn)/^^ Dn is the Fano-factor. 

In this section we show the solution given by our theory for the SRL model, as 
well as how it recovers the the FDT and the NEFDT in the appropriate limits. By 
contrasting these results with the exact solution, we will be able to show that the 
Markovian approximation does not contain quantum fluctuations, thereby needing 
a non- Markovian approach to capture the physics of quantum noise [50]. It is 
therefore interesting to see to what extend the four results coincide, and in what 
regimes our Markovian approach is valid and recovers the proper physics. We will 
see that when kT ^ eV, huj, the theory captures well both the exact solution 
and the FDT and NEFDT. Also, in transport configurations, with the level within 
the bias window, the Markovian approximation agrees well with the exact solution, 
reproducing in particular previous studies with eV ^ kT [50] dSl [HI HZ] • However, in 
a situation with level outside the bias window, the Markovian approach presented here 
does not capture quantum noise physics, effect that we observe at high frequencies 
{huj > kT, eV). Although in this situation transport due to cotunneling processes 
becomes more relevant, the difference is due to the Markovian assumption as can 
be seen using a non-Markovian extension of the theory [SO]. We also study the 
finite- frequency skewness of the current distribution as given by equation (j35|) . This 
shows to be insensitive to thermal fluctuations near equilibrium, therefore revealing 

* For a non- symmetrized version of the noise spectrum through a two-terminal conductor c.f. 1651 . 



+(^--^^oth n^^] 

2 \ 2kT I 
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Figure 1. Zero- frequency noise S^^^ {uj = 0) at finite voltage for the single 
resonant level model as a function of temperature. For comparison we also show 
the FDT, the NEFDT, and the exact solution. In the main figure the level is 
within the bias window £ = 0. The inset shows a regime with the level outside 
the bias window e/eV = 5. Rest of parameters: Vi^/eV = Vn/eV = 0.1. In the 
main figure all the quantities coincide when kT > eV, as expected. In the inset all 
fluctuations are thermal and therefore, all quantities coincide in the whole range of 
temperatures. The typical physical units are T ~ 10 — 100 mK, y ~ 10 — 100 fiV, 
r ~ 10 - 100 MHz. 



the "shot" contribution in a situation in which thermal fluctuations dominate in the 
noise spectrum {kT ^ eV, hw). 

3.1. Zero- frequency counting statistics. 

Let us start by showing the zero-frequency noise spectra corresponding to the SRL 
model. Although this limit has already been studied in detail, a full comparison 
between our theory and the exact solution will help to understand the finite-frequency 
case. Particularly important is the linear response regime, at which studies of this 
model are scarce. As mentioned before, in this regime the noise should exhibit thermal 
fluctuations in order to fulfill a fluctuation-dissipation relation, while the skewness, 
on the other hand, should go to zero as the voltage V goes to zero [55] . Fig. [1] shows 
how our calculation captures correctly the FDT, 5^^^ = 2kTG, in the proximity of 
linear response, kT > eV. For comparison, we also plot the zero-frequency limit of 
the NEFDT in Eq. gT]), namely S'^^) = /^(^jcoth(^) [55]. In the opposite regime, 
kT < eV, the Markovian approximation is larger than the exact solution, discrepancy 
that can be understood as originated from the lack of cotunneling contributions in our 
calculation [7^. As expected, below kT/eV ^ 1 the FDT is not fulfilled. We can also 
see that the NEFDT, exact for tunnel junctions, for a two-terminal device performs 
quite badly when kT < eV, but correctly in the opposite limit. This failure of the 
NEFDT at low temperatures disappears when the level is outside the bias window. 
This is a low-current regime, and thus a tunneling limit. The inset of Fig. [1] shows this 
situation, where all fluctuations are thermal and the four curves coincide in the whole 
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Figure 2. a) Noise and skewncss near linear response, for hVi^/kT = hVn/kT = 
eV/kT = 0.05, as a function of e. b) Noise and skewness as a function of voltage. 
We also show the equilibrium fluctuation-dissipation expression 2kTG and the 
exact solution for comparison. 



range of temperatures. At finite frequencies we will expect a quantum noise step in 
the spectrum at frequencies huj ^ e, effect that will be studied in the next subsection. 

The behaviour of the zero-frequency noise spectrum close to equilibrium with 
respect to e is shown in Fig. [2^. Here we see how the FDT is fulfilled by our theory and 
the good agreement with the exact solution. We also plot the zero-frequency skewncss, 
that although of small magnitude in the same scale, is nonzero in a situation where the 
noise spectrum is completely dominated by thermal fluctuations. This insensitivity 
of the skewness to temperature allows us to extract intrinsic correlation effects in 
near-equilibrium conditions. In Fig. we show the same quantities as a function 
of voltage. Interestingly, the Markovian result coincides with the exact solution in 
the whole range of voltages. The FDT, however, starts to disagree with these for 
voltages eV/e > 0.2. As anticipated, the skewness vanishes as the voltage goes to 
zero. In Fig. we plot noise and skewness as a function of s for increasing voltages. 
As the bias increases, the skewness (dashed lines) shows peaks evolving into plateaus 
at values of e corresponding to the chemical potentials of the reservoirs. This effect, 
which is due to non-equilibrium fluctuations, is completely masked in the noise (solid 
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Figure 3. a) Increasing noise (black full lines going up) and skewness (red 
dashed lines going up) as a function of e for increasing voltages eV/hFji = 
1,10,20,40,60,80,100,120, (parameters kT/hVn = 20, r^/r^ = 1). b) Excess 
noise 5(2) (y) - = 0) and skewness versus e for increasing voltages. 



lines) even at the highest voltages due to thermal fluctuations. This is clearly seen in 
Fig. [31d, where we show the same comparison after substracting thermal fluctuations 
to the noise value (excess noise defined as S<^'^\V) - S'^^^V = 0)). Here it is clear that 
at low detuning, e < eV/2 (position of the peaks in the figure), and when kT > eV, 
the skewness can reveal the "shot" contribution, while this is masked by thermal 
fluctuations in the noise spectrum. 

3.2. Finite-frequency counting statistics. 

To study the case of finite frequencies, we use our formulae ([55|) - ([55)) applied to 
the SRL model. In a situation with the level within the bias window, we find a 
similar behaviour to Fig. [T] However, now the NEFDT - equation (|Tf| - is fulfilled 
for temperatures kT > eV + hut. This is shown in Fig. 01 Remarkably, at finite 
frequencies the Markovian approximation is basically exact in this direct transport 
regime. In the high-bias regime eV ^ huj, kT we also find that the Markovian 
approximation agrees perfectly with the exact result (not shown), in accordance with 
previous studies [80l [151 [HI [17] . When the level lies outside the voltage window, the 
situation changes drastically (see inset of Fig. [J) . Here the Markovian approximation 
is no longer appropriate when kT < eV + huj. Both the exact result and the NEFDT 
contain quantum fiuctuations, while the Markovian calculation only captures the 
thermal contribution (and therefore fulfills the FDT). The exact solution presents 
a small structure at temperatures of the order of e. This cannot be resolved with the 
NEFDT. As expected, when kT > eV + Hlu, all curves coincide. 

The general trends explained so far become more evident in Fig. [SI where we plot 
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Figure 4. Noise, FDT and NEFDT as a function of temperature for u)/eV = 10. 
In the main figure the level is within the bias window e = 0. The inset shows a 
regime with the level outside the bias window e/eV = 5. Rest of parameters: 
hTi^/eV = hFji/eV = 0.1. In the main figure all quantities coincide when 
kT > eV + hoj. In the inset all fluctuations are thermal and, therefore, the 
shot noise and the FDT coincide in the whole range of temperatures. When 
kT < eV + hw, the NEFDT is above due to quantum fluctuations. The 
exact solution contains also corrections due to cotunncling processes, which are 
dominating in this regime. 



noise and the fluctuation-dissipation theorem near hnear response as a function of 
frequency. Here the Markovian noise is always flat and equals 5'^-' = 2kTG, whereas 
the NEFDT and the exact solution he above and show quantum noise steps. Let us 
start by considering the case £ = 0. In the whole range of frequencies, the Markovian 
approximation is basicahy exact in this situation of direct transport. The NEFDT 
shows the correct zero-frequency limit, since then the fluctuations are purely thermal. 
At high frequencies, however, the NEFDT converges to the Poisson value of a single 
barrier with tunneling rate r/2, being F := (F^ -I- Ffl)/2, (c.f. S}^ in the plot). 
This is in agreement with the validity of equation (j47p for a tunnel junction, and in 
contrast with the exact solution, which contains partitioning, and therefore its w — )■ oo 
limit is r/4. In the case in which the energy lies outside the bias window {e/kT = 5 
in the figure), transport is possible because of the finite temperature as well as due 
to quantum fluctuations. The Markovian noise only contains the former and is flat 
with frequency, while the exact result contains both and shows a quantum noise step 
centered at hw/kT = e/kT = 5. Although in this regime cotunneling contributions 
are important, the difference lies in the Markovian approximation, as can be seen 
using a non- Markovian theory [50]. The NEFDT shows in this situation a quantum 
noise step centered at fiu; = 2£. This discrepancy with respect to the exact solution 
can be understood in terms of the tunneling approximation leading to Eq. (j47p . which 
presents a step located at an effective chemical potential 2£. Again, the high frequency 

(2) 

limit coincides with that of S)^ . The intermediate regime where huj ^ eV is studied 
in the inset. Here, we set eV/kT ~ 25 and observe a flat behaviour for the Markovian 
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Figure 5. S'^' near linear response {eV/kT = 0.0005) as a function of 
frequency. For comparison we also show the NEFDT (Eq. II47I I) and the exact 
solution. S^'^'^(u) is flat for the whole range of frequencies, and coincides with 
the equilibrium FDT as expected. The NEFDT however disagrees with those 
two, showing also quantum fluctuations which are absent in the Markovian noise 
spectrum. The quantum noise steps shown by the NEFDT are however at 
hu) = 2£, in contrast to the exact solution, which shows steps at hw = e. This 
is due to the fact that the NEFDT works well for tunnel junctions, but does not 
capture partition noise. This becomes clear also from the saturation value at large 
frequencies, as described in the text. Rest of parameters: /iF^/fcT = hTn/kT = 
0.05. The inset compares the exact solution with the Markovian approximation 
for a different regime, namely eV/kT = 25. We see that while the Markovian limit 
is flat for all frequencies, the exact solution presents a dip a.t twj = it|£ it eV/2\. 
Rest of parameters: e = 0, hVL/kT = hVn/kT = 0.25. 



solution, whereas the exact solution presents a dip at hiu = ±|e± eV/2\ (coinciding 
with the position of the chemical potentials with respect to the energy level). This 
clearly illustrates how the Markovian approximation captures well the physics in the 
linear response regime and a direct transport configuration, but when the frequency 
is comparable to the applied bias, it fails to capture the quantum noise. 

We now proceed to discuss the finite-frequency noise and skewness spectra of the 
total- and particle- current distribution. In the previous discussion, the Markovian 
noise was always flat as a function of frequency, something which is well known 
for symmetric systems (Fl = r_R) - see for instance [5]. In the case Tl 7^ r^, a 
proper partitioning of displacement currents (see discussion in subscction l2.5p becomes 
essential as we will show next, and the way this is made affects significantly the 
spectrum. Fig. [6^ shows S'''^\u}) in a transport configuration, eV/hVj^ ~ 50 and 
e = 0. As in the results shown previously, the total noise spectrum is flat for a 
symmetric configuration. Interestingly, this flat behaviour persists even when the 
system is made asymmetric (Tl 7^ ^b.)- This is due to the current-partitioning model 
assumed here: Itot = chIl + PIr with ex. = Tr/V and j3 = Tl/T. The noise spectrum 
corresponding to particle currents displays information about the rates; in contrast to 
the total noise, it shows a dip with half- width 2T. In Fig. HId we show the skewness 
along the representative direction u' = — w. Interestingly, the skewness corresponding 
to the total current starts to develop a dip that shows the asymmetry of the system. 
The particle-current skewness presents a similar behaviour to the noise counterpart. 
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Figure 6. a) 5(2) (oj) and b) S^^^ui, -oj) for eV/hFji = 50, kT/hVji = 20, and 
£ = 0. The spectra for particle currents and for total currents significantly deviate 
form each other, even for large asymmetry. The insets correspond to noise and 
skewness through the left barrier for F^/Ffl = 10. 



However, for {i- /2 i^T l/T r (3 + \/5)/2 it develops a minimum whose position 
depends on the value of the rates. In the asymmetric case, ^ Tr, the particle- 
current noise and skewness can even present different lineshapes. This can be seen 
contrasting the insets of Fig. [5^. and Fig. ^p. In the linear response regime the curves 
for the noise look similar (not shown). The skewness on the contrary goes to zero 
in magnitude and shows a structure that depends on temperature, and that changes 
from a dip to a peak as e is increased from zero to a finite value. In summary, 
we see that the spectra for total and particle currents differ significantly from each 
other even for large asymmetry. This means that the assumption of calculating noise 
spectra using particle currents only, used commonly in the literature, is flawed. Here 
we have assumed the current partitioning given by a = Fjj/F, /? = Tl/T. If the more 
simplistic partitioning a = /3 = 1/2 is assumed, the results for the total cumulants in 
the asymmetric case change drastically (not shown). In particular, the noise is then 
no longer flat but has a dip structure, and the skewness shows a peak around zero 
frequency. 
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4. Conclusions 

In conclusion, we have developed a theory of frequency-dependent counting statistics 
of electron transport through nanostructures within the framework of Markovian 
quantum master equations. We have illustrated our method with calculations of noise 
and skewness in a single resonant level model at finite bias voltages and frequencies. 
By comparing with both the exact solution and the finite-frequency version of the 
noncquilibrium fluctuation-dissipation theorem, Eq. (j47|) . we have identified the 
regimes of validity of our Markovian theory at flnite frequencies. In particular, wc have 
shown that the Markovian limit is basically exact in transport configurations (level 
within the bias voltage window), as long as Hco ^ eV or fiw <^ eV. In intermediate 
situations where hui ~ eV, or with the level outside the bias window, the Markovian 
limit fails at finite frequencies due to the lack of quantum fluctuations [50] . 

We have also discussed how the noise spectra for particle currents and for total 
currents significantly deviate from each other, even for large asymmetries Tr/^l 1- 
This demonstrates that calculating spectra using particle currents only leads to 
incorrect results in general. Our method allows the calculation of finite-frequency 
current cumulants of arbitrary order, as we have explicitly shown for the second and 
third order cumulants, Eqs. and ([55)1 . These formulae generalize previous zero- 
frequency expressions and can be viewed as an extension of MacDonald's formula 
beyond shot noise. Recently, this has been extended to study the time-averaged shot 
noise spectrum in the presence of periodic ac fields [STJ |55]. Interesting extensions 
of our study along these lines would allow us to study frequency-dependent high- 
order cumulants in nanostructures driven by time-dependent fields, or, even more 
challenging, in systems showing nontrivial non-linear dynamics such as self-sustained 
oscillations without external time-dependent driving |83j . 
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Appendix A. Derivation of frequency-dependent cumulants 

The expressions (|33p - (|35[) follow from derivatives of moment generating functions. 
Performing derivatives of ([32]) we find 




(A.1) 




= (Z1Z2) z^h^^\2J^^^no{zi2)J^''^ 

+ Zi2J^^'>no{z2)noizi2)J^'^ + Jo'"^) 







(A.2) 
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5(?)>(zi,Z2,^3) = {Z1Z2Z3) 9x1^X2^^X3 (^(X3, 2^3)f^(X23, 223)f^(Xl23, 2123] 

= {Z1Z2Z3) z^h^^l 

X {J^^^no{z3)no{z23)J^'^noiz23)flo (Z123 

+ 2J^'^noiz3)noiz23Wzi23)J^^^no{zi23)J^^^ 

+ 4=Z23Jo'^^oiz23)noizi23)J^''^^oizi23)J^^^ 

+ 2z^3'j^'^noiz23)J^^^no{z23)no{zu3)J^'^ 

+ 62:23"^^l23^0^^^o(^123) Jo ^^^0(^123) Jo 

+ j!t'^no{z3)noiz23)no{zi23)Jo^^ 

+ 2z^^'j^'^no{z23)no{zi23)J^^^ 

+ Z^^'J^^^ noiz23)noizi23)J^'^ 

+ 32:23^^123^0^^^0(^123)^0^' 

+ 3Z;73"^2:f23 Jo"^''^o(2l23) Jo^' 



-^23 ^123^0' )' (^-3) 

where Zij := Zi + Zj, Zijk := Zi + Zj + z^. Next we use 

no{z2)no{zi + Z2) = — \no{z2) - no{zi + Z2)] , (a.4) 

Zl I i 
and add the "lesser" part (<) corresponding to negative Laplace frequencies. At this 
point we change to "physical" frequencies uj := iU2 + ■ . ■ + ujn, := wa + . . . + ujn, 
etc., and symmetrize the result. This means adding the expressions corresponding to 
all the possible frequency switchings. In this step wc take into account that 

lim ( + i ) = lim , = 27r(5(a;), (A.5) 

r;->0 \ iuj + rj —luo + rj J Tj^o u)^ + 

where 77 — ?> is a small parameter coming from the "greater" (>) or "lesser" (<) parts. 
Finally, we make use of this energy conservation inherited from the time-translational 
symmetry of the cumulants. We then arrive to the equations (P5)) - p5p used in the 
main text. Importantly, after frequency symmetrization one can realize that the first 
three cumulant formulae are equal to their moment counterparts. 



A.l. Diagrams. 

Interestingly, the results (|A.ip - (|A.3[) given above can be derived following a 
diagrammatic technique, similarly to how this is done with Fcynman diagrams in 
the expansion of the partition function or the S'-matrix. Similarly we write the CGF 
in terms of a series expansion, either in the time domain or in the frequency space. 
To that end we expand each of the ^-dependent propagators in the CGF as a Dyson 
series: 

^ 00 
^ z-yV{x) ^ ^'^'^ ^'^x^oizT ■ (A.6) 

This suggests the use of diagrams of the form given in Fig. lAll (a) . In the frequency 
domain U this rules are: 



tl If we work in the time domain, it is enough to label the propagating lines with the corresponding 
times at the beginning and end of each line (see Fig. lAll (a,)). 



Finite-frequency counting statistics of electron transport: Markovian Theory 



21 



• To each bare propagator rio(5fe) in the expansion we associate a hne with a 

superscript k = X^i^jv+i-fc where N is the order of the cumulant we want 
to obtain. 

• To each jump operator jT^j. in the expansion we associate an encircled cross with 
superscript k. 

The formula for the generating function to a given order can therefore be written 
diagrammatically. For example, to second order we have 

g{x, z) = TsTr{n{xirziMX2, ~Z2)pt'} 

= TsTr{ino{z2) + no{z2)J^Mz2) + . . .) 

X (f7o(zi + Z2) + VLo{zi + Z2)Jx^+xMzl + Z2) + ...) 

(A.7) 

We can then multiply the different terms using diagrams as described above. The 
multiplication of propagators implies joining them together. The result can be 
simplified using the property Jxi+X'i ~ ^x\^xi + + ^xn which diagrammatically 
is denoted as 

12 12 1 2 

(g) = (g)(g) + (g) + (8). (A.8) 

Here, the super-index 12 denotes an associated frequency z\ + Z2 and counting field 
Xi +X2- Next, to arrive to the frequency-dependent moment^ we take derivatives with 
respect to counting fields. Diagrammatically, the derivative d^^, means removing a 
circle with index k. From here we can rewrite the expression analytically. The outcome 
can be simplified using (jA.4[) . and needs to be multiplied by z\Z2 (case TV = 2), coming 
from the Fourier transform of the time derivatives in the frequency domain. We finally 
need to take the average in the stationary state and symmetrize the result as dictated 
by 7s. 

With the diagrammatic approach we realize that the diagrams contributing to 
the final result can be arranged in tables (see Fig. lAll (b) to (d)). These reproduce 
the results given in (|A.l|) - (|A.3p . To construct these tables we proceed according to 
the following rules: 

• To arrive to an expression for the cumulant of order iV, write a table with iV time 
(frequency) intervals and corresponding superscripts k. The propagation of time 
will be taken from right to left. 

• Write all the possible diagrams having crosses (jumps) distributed in the 
different intervals, with the constraint that the maximum number of crosses in 
each is set by the corresponding index k. Diagrams with n jumps occurring at 
the same time have to be included as well. These crosses are enclosed together 
with a box, and contribute with the jump operator := d^J-^^^Q. 

• Taking into account that jumps occurring in the same interval are indistinguish- 
able, and that each cross can be associated to one of the possible counting fields 
Xi, . . . , Xfe present in that interval, write the multiplicity of each diagram on the 
right. 

• Write the mathematical expression corresponding to each diagram (see Fig. lAll 
(a)) and sum the different terms evaluated at z = icj. 

• Take the average in the stationary state, and multiply by {—i)^. The resulting 
expression corresponds to the unsymmetrized ("greater", >) moment of the 
number of particles. 
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• Multiply by (iioi) . . . (iwjv). This gives the unsymmetrized moment of the current 
distribution. 

• Add the "lesser" (<) part, that is, the expression corresponding to negative 
frequency. 

• Finally, symmetrize the result, adding all the possible switchings of frequencies. 
This gives the symmetrized moment of the current distribution. The result can 
be simplified using (jA.4|) and ()A.5j) . 

As mentioned above, explicit derivation gives the same result for the expressions 
of cumulants of the current distribution as those derived for the moments up to TV = 3. 
To higher orders it is unknown for us if this property still holds or not. Expressions 

(2) 

for the cross correlations, e.g. S)^^ := (/(ti)/(t2))c, between two (or more) stochastic 
processes, e.g. L and i?, can also be derived with this technique. To this end we 
simply need to label each of the jumps occurring at L or i? accordingly (see Fig. lAll 
(c)), having two types of jump operators, J7l and jT/j. Also, expressions for the total 
current {uIl + PIr) or accumulated current [II — Ir) can be derived using the jump 
operators p9 |) -(|40 |) in the diagrams. 

A. 2. Equivalent form. 

We can write down an equivalent form to expressions (j33])-(|35]). This will allow us to 
obtain an analytical expression for their zero- frequency limit, which is not well defined 
in the form given above. To this end we make use of the projectors P := |0))((0| 
and Q := 1 — P, where P projects onto the subspace spanned by the stationary 
state |0)) = pg*°* lUJ; and we define the pseudo-inverse Ro{z) := Qno{z)Q, such that 
rio(z) = Ro{z) + P/z. Making this change in (jA.l[) - (|A.3p and symmetrizing the 
expression (including positive and negative frequencies) we get 

z(/(z))=<5(z)(J-«), (A.9) 

(A.IO) 

, Z2, Z3) = S{zi + Z2 + Z3) 

X (Jo^'^ + J^'^ [Ro{zi) + Rn{z2) + Ro{z3)] 

+ Jo''^ [Ro{zi2) + i?o(z23) + R^iziz)] Jo^'' 

+ Jo^'^i?o(zi) J-J'^ [Ro{zi2) + Roizis)] Jo^'^ 

+ Jo''^^o(^2) Jo'" [^0(^12) + i?o(^23)] Jo^'^ 

+ J^'^ R0{Z3)J^'^ [Ro{zi3) + Ra{z23)] Jo^'^ 

+ Z-\J^'^)J^'^ [i?0(^12) 

--Ro{z2)+Roizi3)-Roiz3)]J^''^ 

+ Z2\J^'^)J^'^ [i?0(zi2) 

-Roiz,) + Ro{z23) - Ro{z3)]J^'^ 

ft The state ((0| denotes the left eigenvector of the Liouvillian. The tilde indicates that it is not the 
adjoint to |0)), since the Liouvillian is not Hermitian. 
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-Ro{zi) + Ro{z2z) - Ra{z2)] Jo^'^>. (A.ll) 

Now we make use of the delta function to write Z2 = — zi in the noise expression and 
Z3 = — zi — Z2 in the skewness resuh. Performing the change of variables zi — > —iuj, 
Z2 — >■ iio in the noise and zi —ico, Z2 ^ ioj — iuj' , Z3 — >■ iuj' in the skewness, we 
obtain 

ilstat = iJ^'), (A.12) 
I'S^^Hu) = (jj^) + J^''>Ro{iLo)J^'^ + J^'^Roi^iLo')J^'^), (A.13) 
w') = (J-f ) + J-Ji' [Roi-tio) + RoiiLo - tij') + RoitLo')] Jo*'^ 

+ J"J''i?o(-*^) Jo''' [i?o(-*'^') + Roiiu}' - Jj'^ 

+ J-J''i?o(jw - Jo^'' [i?o(-«^') + Rniioj)] Jo''' 

+ Jc|''i?o(jw') Jj'' [^o(«^^' - luj) + Roiioj)] jj'' 

+ {~iLu)-\j^'^)J^'^ [i?oH^') 

-i?o(iw - iw') + Ro{iuj' - iw) - Ra{iuj')] J^^^ 

+ {iLO - zc^')"'(jj"> Jo*'' [R^i-iu:') 

-Ro{-iuj) + i?o(jw) - Ro{i^')] Jj'' 

+ (jc^')"'(jj">jj" [i?oK-«c^) 

-Ro{-ioj) + i?o(jw) - i?o(j'^ - ii^')] Jj'')- (A.14) 

The limit w — !■ of these expressions is well defined, and they can therefore be used 
to check that the proper result is recovered in that limit. 

A. 3. Zero-frequency limit. 

As mentioned, expressions (jA.12[) - (jA.14p are well behaved when a; — > 0. The zero- 
frequency noise comes straightforwardly from (jA.13[) setting a; = 0. For the skewness, 
this limit requires nevertheless noticing that 

lim [Ro[ioj) - i?o(-«w)] = 2ia;9aj-Ro(*'^)U=o- (A. 15) 

oj— >0 

So the zero-frequency skewness can be written as 

z5(3)(0,0) = (jj" +3jJ''i?o(0)Jo''' +3Jo'''i?o(0)jJ" 
+ 6Jo'''i?o(0)ji"i?o(0)Jo'") 

+ 6(Jo''')(J-J''a„i?o(0)jJ"). (A.16) 
Now, since 9^i?o(0) = i?(0)i?(0), we have 

z5(3)(0,0) ^ (Jo'"' +3Jo'''i?o(0)jJ" +3Jo*"i?o(0)Jo'" 
+ 6Jo(')i?o(0)Jo'''i?o(0)Jo'") 

+ 6(J-J'')(jJ"i?o(0)i?o(0)jJ"). (A.17) 
Which is the zero- frequency limit found in |61) . 
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Figure Al. Diagrammatics to obtain the frequency-dependent cumulants. 

a) Building pieces for the diagrammatic technique. A hne is associated to 
a bare propagator, a cross with a jump operator, and a circle with the 
time dependence of J7 (term (e*^ — 1) in the single-particle unidirectional 
tunneling case). Derivatives with respect to counting fields eliminate 
circles correspondingly. Diagrams can be simplified using rules like IIA.81 1. 

b) Diagrams for the noise. Reading this table we find the expression 



(2) (2) 
fto{zi)Qo(z2),jQ iio(-22). c) Diagrams for the second order cross correlation S}^^. 

Reading this table we find the expression fto{zi)J^Q^^^o{zi)^o{z2):jQ^^^oiz2) + 



r(l) 



na{zi)no{z2)j^^lno{z2)j^llno{z2) + no{zi)^o{z2)J^%^o{z2)j^'lno{z2)- d) 

Diagrams to derive the frequency-dependent skewness formula. 
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Appendix B. Derivation of the self-energy 

To calculate the kernel of equation (|T7| . we follow the perturbative treatment by 
Schoeller and coworkers [68l|69]. Let £s, and Ct be the corresponding Liouvillians 
to ©-(HI). The last can be written in the form 

Ct^-i V.»„G«?4S(x), (B.l) 

where f = + (— ) refers to the creation (annihilation) of particles in the leads, and 
G^^ and J|a(x) ^^'^ system and reservoir superoperators respectively, that act on an 
operator A as 

I P = + (B.2) 

Where g^^ and are defined as 

3+ =^(a|d„Ja')|a)(a'|, (B.4) 

aa' 

3Ux)=<^W''''-'\ (B.5) 

and g.^j = {g^)\ j^aix) = (jr^alx))^- The index Sa = ±1 is taken according to the 
sign convention for the current flow in lead a. 

With this notation, the self energy to order |Vp reads [551 [Ml [70] : 

I](2)(z,x)P5(io) 

^ ^TrR,{v,.„G^„^J,^g(x) ^ _ ^^ V,'.'n.'GiO^;tix}p{to)} 



— V /" r^^^' ,(e y)G^p 



-P' 



\a)){{a\G~yfi^Cpe/kT)depsito) 



z - i^{e + Ha) - iK 

[1 ip 
b(AQ + CMa - iz)) + —(j) ip{Xa + i^-a - iz)) J Ps(^o), 

(B.6) 

where the summations run over all scripts, and _D is a high-energy cutoff set by the 
bandwidth of the Fermi leads - larger than the rest of energy scales in the problem. 
In this expression, we have introduced a complete set of eigenstates of the system 
Liouvillian, /3s|i)) ~ iAa|a)), and the definitions f{x) {^e^/^^ + and 

^'Zm'i^^x) ■■= r_„,(£)e^-«("^)'^", (B.7) 

being Tamm'ie) = x Sr; ^rjam^vam'Sie - Er^a) (wlfich wc take to bc independent of 
the energy Tamm'i^) ~ ^amm'), and \1/ is the digamma function. The self-energy 
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(|B.6p is important as it permits us to explore correctly the low bias limit {eV < kT) 
to sequential tunneling order. This self-energy is non-Markovian as the Markovian 
approximation has not been made up to this point. This can be made (together with 
the secular approximation) taking the limit z — > of (|B.6|) . and it is what has been 
used throughout the text. 
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